/*******************************************************************************
																				
	DESCRIPTION:  	This do file produces graphs that depict the heterogeneity 
					by cyclicality in our sample.

*******************************************************************************/

clear all
global id_code 114_5
pause on
set seed 2110


* Import cyclicality data from R
use "${data}/114_3_CyclicalityBetas_Full_2006.dta", clear

* Generate quantiles in terms of alfa and beta from the regression in logs
xtile quintBetaU = beta_u, nq(5)
xtile quintBeta0 = beta_0, nq(5)

xtile decBetaU = beta_u, nq(10)
xtile decBeta0 = beta_0, nq(10)

* Generate empty variables to store the results:
gen u_t = .

foreach i in 0 U {
	forval j = 1/5 {
		gen f_q`i'_`j' = .
	}
	
	forval j = 1/10 {
		gen f_d`i'_`j' = .
	}
}

gen exp = .
gen f_ave = .

* Loop over which coefficients to use:
foreach coef in _shrunk_ind {
	
	* Calculate E(exp(beta_0)) for the left panel:
	replace exp = exp(beta_0`coef')
	sum exp if !missing(beta_u`coef')
	local mean_exp_beta_0 = r(mean)

	* Store 2006 unemployment rate as local:
	local u_2006 = shareUnempSS[1] 

	* Calculate E(exp(beta_0)) x E(exp(beta_u x (u_t - u_2006)) | Q) using 
	* quintiles of beta_u for the LHS panel:
	local i = 1

	foreach u_t of numlist `u_2006' 4.4(.2)8 {
			
		replace exp = `mean_exp_beta_0' * (`u_t' / `u_2006')^(beta_u`coef')
		
		forval q = 1/5 {
			sum exp if quintBetaU == `q'
			replace f_qU_`q' = r(mean) in `i'
		}
		
		forval d = 1/10 {
			sum exp if decBetaU == `d'
			replace f_dU_`d' = r(mean) in `i'
		}
		
		replace u_t = `u_t' in `i++'
	}

	* Calculate E(exp(beta_0 + beta_u x (u_t - u_2006)) | Q) using quintiles of beta_0 for
	* the RHS panel:	
	local i = 1

	foreach u_t of numlist `u_2006' 4.4(.2)8 {
			
		replace exp = exp(beta_0`coef') * (`u_t' / `u_2006')^(beta_u`coef')
		
		forval q = 1/5 {
			sum exp if quintBeta0 == `q'
			replace f_q0_`q' = r(mean) in `i'
		}
		
		forval d = 1/10 {
			sum exp if decBeta0 == `d'
			replace f_d0_`d' = r(mean) in `i'
		}
		
		replace u_t = `u_t' in `i++'
	}

	* Finally, produce overall average profile E(exp(beta_0 + beta_u x (u_t - u_2006))):
	local i = 1
	foreach u_t of numlist `u_2006' 4.4(.2)8 {
			
		replace exp = exp(beta_0`coef') * (`u_t' / `u_2006')^(beta_u`coef')
		sum exp
		replace f_ave = r(mean) in `i'
		replace u_t = `u_t' in `i++'
	}
	
	* Put averages in new frame:
	frame put u_t f_q* f_d* f_ave, into(cycl`coef')
	frame change cycl`coef'
	keep if u_t != .
	
	* LHS plot:
	graph twoway (line f_qU_1 u_t, color(orange_red*0.4) lwidth(medthick)) ///
		(line f_qU_2 u_t, color(orange_red*0.7) lwidth(medthick)) ///
		(line f_qU_3 u_t, color(orange_red*1) lwidth(medthick)) ///
		(line f_qU_4 u_t, color(orange_red*1.3) lwidth(medthick)) ///
		(line f_qU_5 u_t, color(orange_red*1.6) lwidth(medthick)), ///
		yline(`mean_exp_beta_0', lpattern(dash) lcolor(black)) ///
		xtitle("Unemployment rate (%)") xscale(range(4 8)) xlabel(4(1)8) xtick(4(0.5)8) ///
		legend(order(1 "Q1" 2 "Q2" 3 "Q3" 4 "Q4" 5 "Q5") ///
			position(6) rows(1) ///
			symxsize(*0.5) size(medsmall) linegap(2) ///
			region(margin(small)) subtitle(" ", size(minuscule)) ///
			title("Quintiles of {&beta}{sub:U}, Fixing {&beta}{sub:0}", ///
				size(medsmall) alignment(top) color(black))) ///
		ytitle("") ylabel(0.2(0.1)1) ///
		/// t2title("{bf: E}[{&pi}{sub:u}|{bf: E}({&beta}{sub:0}), Q{sub:U}]") ///	
		/// t2title("A. By Quintiles of Cyclicality, Fixing Initial JF") ///
		graphregion(color(white)) ///
		name(g1, replace)

	* RHS plot:
	graph twoway (line f_q0_1 u_t, color(ebblue*0.4) lwidth(medthick)) ///
		(line f_q0_2 u_t, color(ebblue*0.7) lwidth(medthick)) ///
		(line f_q0_3 u_t, color(ebblue*1) lwidth(medthick)) ///
		(line f_q0_4 u_t, color(ebblue*1.3) lwidth(medthick)) ///
		(line f_q0_5 u_t, color(ebblue*1.6) lwidth(medthick)), ///
		yline(`mean_exp_beta_0', lpattern(dash) lcolor(black)) ///
		xtitle("Unemployment rate (%)") xscale(range(4 8)) xlabel(4(1)8) xtick(4(0.5)8) ///
		legend(order(1 "Q1" 2 "Q2" 3 "Q3" 4 "Q4" 5 "Q5") ///
			position(6) rows(1) ///
			symxsize(*0.5) size(medsmall) linegap(2) ///
			region(margin(small)) subtitle(" ", size(minuscule)) ///
			title("Quintiles of {&beta}{sub:0}", ///
				size(medsmall) alignment(top) color(black))) ///		
		ytitle("") ylabel(0.2(0.1)1) ///
		/// t2title("{bf: E}({&pi}{sub:u}|Q{sub:0})") ///
		/// t2title("B. By Quintiles of Initial JF") ///
		graphregion(color(white)) ///
		name(g2, replace)
		
	* Combine graphs:
	graph combine g1 g2, rows(1) ycommon ///
		l1title("Average 6-Month JFR", size(medium)) ///
		graphregion(color(white)) ysize(4.5) xsize(8) name(g3`coef', replace)
	
	* Store them:
	if "`coef'" == "_shrunk_pop" {
		graph export "${output}/${id_code}_Visualizing_Heterogeneity_Cyclicality_ShrunkPop.pdf", as(pdf) replace 
	}
	else if "`coef'" == "_shrunk_ind" {
		graph export "${output}/${id_code}_Visualizing_Heterogeneity_Cyclicality_ShrunkInd.pdf", as(pdf) replace 
	}
	else {
		graph export "${output}/${id_code}_Visualizing_Heterogeneity_Cyclicality.pdf", as(pdf) replace 
	}
	
	frame change default
}

* Plotting 1-JFR instead (only RHS plot, with deciles, for shrunken coefficients):
frames change cycl_shrunk_ind 

forval i = 1/10 {
	gen pi_d0_`i' = 1 - f_d0_`i'
}

forval i = 1/5 {
	gen pi_q0_`i' = 1 - f_q0_`i'
}

frames change default


* Plotting magnitudes relative to average (only RHS plot, with quintiles, for 
* shrunken coefficients):
frames change cycl_shrunk_ind

forval i = 1/5 {
	gen rel_f_q0_`i' = f_q0_`i' / f_ave
	gen rel_pi_q0_`i' = pi_q0_`i' / (1 - f_ave)
}

graph twoway (line rel_f_q0_1 u_t, color(ebblue*0.4) lwidth(medthick)) ///
	(line rel_f_q0_2 u_t, color(ebblue*0.7) lwidth(medthick)) ///
	(line rel_f_q0_3 u_t, color(ebblue*1) lwidth(medthick)) ///
	(line rel_f_q0_4 u_t, color(ebblue*1.3) lwidth(medthick)) ///
	(line rel_f_q0_5 u_t, color(ebblue*1.6) lwidth(medthick)), ///
	yline(1, lpattern(dash) lcolor(black)) ///
	xtitle("Unemployment rate (%)") xscale(range(4 8)) xlabel(4(1)8) xtick(4(0.5)8) ///
	legend(order(1 "Q1" 2 "Q2" 3 "Q3" 4 "Q4" 5 "Q5") ///
		position(6) rows(1) ///
		symxsize(*0.5) size(small) linegap(2) ///
		region(margin(small)) subtitle(" ", size(minuscule)) ///
		title("Quintiles of {&beta}{sub:0}", ///
			size(small) alignment(top) color(black))) ///		
	ytitle("Relative 6-Month JFR (Average = 1)") ylabel(0.5(0.25)1.75) yscale(range(0.4 1.8) titlegap(5)) ///
	/// t2title("{bf: E}({&pi}{sub:u}|Q{sub:0}) / {bf: E}({&pi}{sub:u})") ///
	graphregion(color(white)) ///
	name(rel_f, replace)


graph twoway (line rel_pi_q0_1 u_t, color(ebblue*0.4) lwidth(medthick)) ///
	(line rel_pi_q0_2 u_t, color(ebblue*0.7) lwidth(medthick)) ///
	(line rel_pi_q0_3 u_t, color(ebblue*1) lwidth(medthick)) ///
	(line rel_pi_q0_4 u_t, color(ebblue*1.3) lwidth(medthick)) ///
	(line rel_pi_q0_5 u_t, color(ebblue*1.6) lwidth(medthick)), ///
	yline(1, lpattern(dash) lcolor(black)) ///
	xtitle("Unemployment rate (%)") xscale(range(4 8)) xlabel(4(1)8) xtick(4(0.5)8) ///
	legend(order(1 "Q1" 2 "Q2" 3 "Q3" 4 "Q4" 5 "Q5") ///
		position(6) rows(1) ///
		symxsize(*0.5) size(medsmall) linegap(2) ///
		region(margin(small)) subtitle(" ", size(minuscule)) ///
		title("Quintiles of {&beta}{sub:0}", ///
				size(medsmall) alignment(top) color(black))) ///	
	ytitle("Relative LTU Risk (Average = 1)") ylabel(0.5(0.25)1.75) yscale(range(0.4 1.8) titlegap(5)) ///
	/// t2title(" [1 - {bf: E}({&pi}{sub:u}|Q{sub:0})] / [1 - {bf: E}({&pi}{sub:u})]") ///
	graphregion(color(white)) ///
	name(rel_pi, replace)
	
grc1leg2 rel_f rel_pi, rows(1) legendfrom(rel_f) ///
	graphregion(color(white)) ysize(4.5) xsize(8) name(rel_all, replace)

graph display, ysize(4.5) xsize(8)

graph export "${output}/${id_code}_Visualizing_Heterogeneity_Cyclicality_ShrunkInd_Rel.pdf", as(pdf) replace 

frames change default
	
* Store desired variables into frame:
frames put LopNr_PersonNr InLnr p_emplAft6M_0M_In beta_0 beta_0_shrunk_ind quintBeta0 ///
	beta_u beta_u_shrunk_ind quintBetaU, into(data2)
	